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Abstract 

The concept of Self-Organized Criticality (SOC) was proposed in an attempt to explain 
the widespread appearance of power-law in nature. It describes a mechanism in which 
a system reaches spontaneously a state where the characteristic events (avalanches) are 
distributed according to a power law. We present a dynamical systems approach to 
Self-Organized Criticality where the dynamics is described either in terms of Iterated 
Function Systems, or as a piecewise hyperbolic dynamical system of skew-product type. 
Some results linking the structure of the attractor and some characteristic properties of 
avalanches are discussed. 

1 Introduction 

The notion of Self- Organized Criticality (SOC) has become a new paradigm to explain the 
widespread appearance of power-law in a multitude of examples like the distribution of the 
size of earthquakes, 1/f- noise, amplitudes of solar flares, species extinction .... to name only a 
few cases pfl E] ■ In this paradigm, the dynamics occurs as chain reactions or avalanches. A 
stationary regime is reached where the distribution of avalanches follows a power law, namely 
there is scale invariance reminiscent of thermodynamic systems at the critical point. A local 
perturbation can induce effects at any scale and there are long-range spatial and time cor- 
relations. In other words, in this paradigm the system reaches spontaneously a critical state 
without any fine tuning of some control parameter. Several models have been proposed like 
the sandpile model 0IH], the abelian sandpile [TO] or the continuous energy model [Hj. The 
results available are mainly numerical and only a few rigorous results are known. Moreover, 
the scientific community is still lacking a general formalism to handle these models properly. 
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It is natural to try to access the macroscopic behavior of large sized systems from the micro- 
scopic dynamical evolution. In this spirit, we have developed a dynamical system description 
for a certain class of SOC models (like Zhang model |14j). for which the whole SOC dynamics 
can either be described in terms of Iterated Function Systems, or as a piecewise hyperbolic 
dynamical system of skew-product type. Several deep results from the theory of hyperbolic 
dynamical systems can then be used, having interesting implications on the SOC dynamics. 
The general setting and some of these results are discussed in this paper. 

2 Definitions. 

In this paper we consider the Zhang model which is defined as follows. Let A be a d-dimensional 
box in 7L d , taken as a cube of edge length L for simplicity and let 9 A be the boundary of A, 
namely the set of points in 7L d \ A at distance 1 from A. Call N = $A = L d , where denotes 
the cardinality of a set. Each site % G A U 9 A is characterized by its "energy" X;, which is a 
non-negative real and finite number. An energy configuration is a vector X = {Aj}. gA . The 
sites of 9 A have always zero energy : as discussed below this mimics energy dissipation at the 
boundaries. Let E c be a real, positive number, called the critical energy, and M. = [0, E C [ N . An 
energy configuration X is called stable when X G M. and unstable otherwise. In an unstable 
configuration the sites i, such that Xi > E c , are called active or unstable. The dynamics on X 
depends whether X is stable or unstable. 

If X is stable, one chooses a site % G A at random with probability ■h, and add to it the energy 
5=1 (excitation). If X is unstable, each active site loses a part of its energy, redistributed 
in equal parts to its 2d neighbors in the following way (relaxation). Fix e G [0, 1[ and set 
a = When i is active it gives the energy aXi to its 2d neighbors and keeps the energy 

eXj. Therefore, the energy is locally conserved during relaxation. The relaxation dynamics is 
synchronous and it is useful to express it in terms of the map: 

F(X) = X + aA [Z(X) * X] . (1) 

where Z(X) is a N dimensional vector, such that -Zi(X) = if Xi < E c and Zj(X) = 1 if 
Xi > E c . The * denotes the product component by component: if X, Y are N dimensional 
vectors, X * Y is the N dimensional vector of components X^Yi. A is the discrete Laplacian on 
A with zero boundaries conditions. Due to the presence of thresholds, F corresponds to energy 
propagation by singular diffusion. 

Zhang model dynamics consists therefore in resting periods where energy is injected and 
stored locally (excitation), followed by periods of activity called avalanches: when the energy 
of a given site exceeds the threshold, this energy is partially released and distributed to the 
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neighbors. This leads to a chain reaction that can propagate on very long scales. However, the 
energy reaching the boundaries is lost since we imposed the boundaries to have always a zero 
energy. This mimics a dissipation mechanism which implies that all avalanches stop within a 
finite number of iterations. Note that with these updating rule, each avalanche starts from only 
one active site. 

The structure of an avalanche is encoded by the sequence of active sites: 

c = {Ct}i< t < TC , 

where: 

Ct = {j G A\Xj > E c in the t-th step of avalanche} , 

and where tq, called the avalanche duration, is the smallest positive integer such that C TC +i = 0- 
Each avalanche can be labeled by a double index (i,j). The first index refers to the site 
where the energy is dropped and the second index labels the different avalanches starting at i 
(including the "empty" avalanche where the excitation of % does not render it active). Moreover, 
to each avalanche corresponds a convex domain M.Uj) in A4. For each i, the domains 
M-uj) form a partition of (HJ- 

The total energy of a stable configuration in a finite lattice being finite, (it is bounded by 
L d E c ), the total number of different avalanches is finite for finite L, (but diverges as L — » oo). 
The size s of an avalanche is the total number of active sites, s(C) = #UtCt- The area a is the 
number of distinct active sites in C. We will generically denote by n an avalanche observable 
(size, duration, area) and n(i,j) will be the value that n takes in the avalanche 

It is numerically observed that, after a sufficiently long time, the dynamics reaches an out of 
equilibrium stationary state 1 where the energy injected per unit time is equal, on average, to the 
energy dissipated at the boundaries. In this regime, the avalanche observables are distributed 
according to a truncated power law: 

P L (n) = ^h(n), 1 < n < & (2) 

Here and in the sequel, the subscript L will refer to the size of A. Cl is a normalization 
constant depending on L, and £l is the maximal value that the observable can take in a box 
of size L (note that this quantity depends also on E c ,e). As discussed above ^ is finite for 
finite L, but it diverges (typically like L@) when L — > oo. r is called the critical exponent of 
the avalanche observable, /z(n) is a finite size cut off, accounting for boundaries effect, and 
such that lim^oo /z,(n) = 1, Vn < Consequently, as L — > oo, Pl(ti) converges to a power 
law distribution P*(n) = 4r. 

1 A precise definition will be given below. 



3 



3 Avalanche maps. 



Each avalanche maps a stable configuration X G -Muj) to the next stable configuration 
obtained when injecting the energy 5 at site i. Consequently, to each avalanche we 
associate a map : A4(i,j) — ► .M such that : 

T M (X) = F^X + fei), X g M(ij), (3) 

where is the z-th canonical basis vector of R^. Note that Tqa is a simple translation in the 
case where the excitation of i does not render it active. 



Let A be the set of all avalanches symbols The transition — * (12, 32) is fega/ iff: 

T (h,h) [M( ildl )]nM iiad2) ?Q (4) 

This means that some energy configuration X G jW^jj) can undergo the avalanche (12, jz) 
after the avalanche (ii, ji) . Note that since the excited sites are selected randomly and in- 
dependently, there is no constraint in the choice of the symbol %2 ■ In other words, there are 
at least iV legal transitions for each symbol The avalanche transition graph Qj± is the 

graph of legal transitions. It is tempting to use the avalanche coding for symbolic coding of the 
dynamics. However, the coding is in general not one to one. 

The maps Tuj\ have some remarkable properties 0: 

(i) Each Tuj\ is affine: 

T (m)( X ) = %,j)( X ) + C (.i,J) ( 5 ) 
where Luj) is an N x N matrix and Cu^ a constant. 

(ii) Each Tuj\ is a quasi- contraction: 

p(L { i,j)) < 1 (6) 

where p is the spectral radius. Moreover, the number of eigenvalues of modulus strictly 
lower than one is equal to the avalanche area a(i,j) 

(iii) The determinant of Luj\ (local volume contraction) is directly related to the avalanche 
size s(i,j) by: 

detL (lJ) = e s ^ (7) 
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(iv) There exists t = t (L, E c , e) < oo and r] < 1, such that, for any legal path j tQ of length 
to on Qj± one has: 

^h ])<^< 1 (8) 
where i is the composition of the matrices along 7^. Namely, each composed 
map TV 1 is a contraction. 

It follows from these properties that Zhang model can be described as a probabilistic graph 
directed iterated function systems. 



4 Dynamical system. 

Zhang model can also be described as an hyperbolic dynamical system of skew product type, 
with singularities. The phase space is the set Ai = E^ X Ai, where E^ is the set of right 
infinite excitation sequences a = {ai, . . . , a*,, . . . \dk € A} such that a t is the £-th excited site. 
A point in Ai is denoted by X = (a, X). We call tt + (resp. 7r~) the projection on E^ (resp. on 
M) such that 7r + (X) = a and tt~(X) = X. 

The dynamic evolution is given by the map T : M. — » M. such that: 

T(X)^(aa,T 8l (X)) (9) 

where a is the left shift on H\. T ai is the mapping from Ai — > whose restriction to Ai ai ,k 
is 7( aiife ). 

The shift a is conjugated to Nx mod 1 on [0, 1]. Therefore, the dynamical system Q has 
a positive Lyapunov exponent Al(0) = log(N). Due to the property (jHj) it has also strictly 
negative Lyapunov exponents Xi{i), i = 1 . . . N corresponding to the projection of the dy- 
namics on Ai. The maps T ai are discontinuous at the boundaries of the domains Ai ai ,k- 
The union of these boundaries, S = U ai=1 iv uj^ 1 ^ dAi ai ,k is called the singularity set. This 
is the set of energy configurations such that at least one site has an energy exactly equal to 
E c . Therefore, the slightest change in its energy can change dramatically the further evolution. 
The singularity set plays therefore an important role (see [S] for a discussion). 

As discussed above Zhang model is expected to reach a stationary state when time tends 
to infinity. In the context of the dynamical system Q it is natural to define this state as the 
Sinai-Ruelle-Bowen (SRB) measure. This is the weak limit: 

fa = lim T*'(/i) (10) 

t— >oo 
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where /i is the Lebesgue measure on M.. Note that the existence of the limit (jl(J|) in Zhang 
model can be rigorously established only in some restricted case jH]. However this is a natural 
assumption from a physicist point of view, since it means that the stationary regime does not 
depend on the way the system is prepared. It is supported by numerical simulations and also 
by some rigorous results established in |U |S] . 

5 Symbolic coding. 

As discussed above, the avalanche coding is not suitable for symbolic coding. However, it is 
possible to construct a finer partition allowing a one to one correspondence between an infinite 
sequence of symbols and a point in A4. More precisely, it can be shown, in some cases, that 
the dynamical system Q has a finite Markov partition. We conjecture 2 that this is true for 
generic values of E c . 

Call V = \V U \ this partition, where f2 is the finite set of symbols parameterizing 

the Markov partition elements.. By construction each element has a product structure 
7r + (Vuj) x 7r~(Vu)- The projection tt + (Vu}) of V u on is a cylinder set consisting of all 
sequences starting with a definite a\ G A, while the projection tv~(^PJ) C M.. For simplicity, 
we will use the notation V w = Tr~(Vu,) 

To each element u corresponds a sub-domain V w C Aiij and henceforth a unique avalanche 
(z,j). Equivalently, to each u one can associate a double index (i (uj) ,j (a;)) where the first in- 
dex refers to the site where the avalanche starts and the second to the corresponding avalanche. 
Note however that, in general, several symbols correspond to the same avalanche since a domain 
Aiij is usually composed by several partition-elements V w . However, in order to simplify the 
notation we will identify uj and the avalanche (i(u) , j (uj)) whenever this causes no confusion. 
To summarize, we have now a symboling coding useful both for the dynamics and for labeling 
the avalanches. 

Moreover, the evolution of a probability distribution on Q can be encoded in Markov transition 
graph with a transition matrix W such that: 

(11) 

2 This conjecture is based on the following (numerical) observation [S]. The £il measure of the ^-neighborhood 
of the singularity set S decreases like u a where a < 1. From this property and from the Borel-Cantelli lemma 
it follows that Lebesgue almost every point in M. has a local stable manifold of positive diameter. Each point 
has also an unstable manifold (corresponding to the direction of the shift in the extended phase space). This 
implies the existence of a finite Markov partition. 
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It can be shown in some cases that W is irreducible jH] . As discussed above (see the comments 
about the weak limit (|1U|)) it is natural to expect that this is a generic situation. Consequently, 
the corresponding Markov process admits a unique invariant measure uil such that tul = itilW. 
The induced measure on the space of legal right infinite sequence on Q, Eq, is conjugated to 
the SRB measure (HUJ) 0. 

6 Main results. 

We have now two possible pictures for Zhang model. This is a probabilistic graph directed 
iterated function system with invariant probability m^, where each map T w has a volume 
contraction e s ^\ Note however that the maps are not conformal. Under the assumption of 
irreducibility of the graph, there is a unique attractor. As discussed below it is expected to be 
fractal for sufficiently small E c values. Zhang is also an hyperbolic dynamical system of skew 
product type, with singularities. Both aspects are fruitful. In this section, we discuss the main 
results obtained from the formalism developed above. 

6.1 Average contraction rate and avalanche size. 

It follows from eq. (J7|) and the ergodic theorem that : 

N 

EM0 = l°g(e)<«> £ - (12) 
i=i 

where (s) L is the average avalanche size. There is therefore a strong connexion between the 
microscopic volume contraction and the average avalanche size. 

More generally, it follows from eq. (JJJ) and (J2J) that the local volume contraction of the maps 
T w on a typical trajectory is distributed according to a truncated power law of type (J2J). This 
suggests therefore a strong connexion between the fractal properties of the stationary state and 
the avalanche size distribution (jSlEj)- 

6.2 Hausdorff dimension versus E c . 

The following proposition can easily be proved. 

Proposition 1 fx^ is singular for all E c sufficiently small. 

Indeed, the L\ norm of all map is bounded from above by one minus the energy flux 
dissipated at the boundaries [7j. This flux e [0, 1] and tends to 1 when E c — ► 0. Consequently, 
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for E c << 1 one can make the Li norm of all map T u arbitrary small. Since the expansion 
is constant it follows that det(DT) < 1 hence p,L is singular (but, as a SRB measure, it is 
absolutely continuous along the unstable foliation). 

When E c increases, the average contraction rate decreases since the smaller is the average 
avalanche size < s >l- Therefore, for fixed N, there exists an E*(N) which is the unique E c 
value such that: 

log(e) < 8 > L + log(iV) = (13) 

For E c < E*(N) the contraction dominates the expansion, while it is the opposite for E c > 
E*(N). Clearly, the invariant set structure is different in these two cases. On the one hand, 
for small E c values, the images of the domains are thin bands which are stretched slower 
than they contract. Therefore, the invariant set has a Cantor structure with large gaps. On the 
other hand, when E c > E*(N), the successive images of the domains P w fill more and more the 
phase space. Consequently, the Hausdorff dimension of the attractor is expected to increase for 
increasing E c , E c < E*(N) and is likely to be constant when E c > E*(N) . 
We conjecture the following: 

Conjecture 1 The Hausdorff dimension of \xl is piecewise continuous and monotonously in- 
creasing on the domains of continuity. 

This is supported by the following argument. On open intervals Jj of E c the structure of 
the mappings T^k) does not change but the domains of continuity M.{i,k) do. Furthermore for 
E c decreasing the probabilities for avalanches with higher contraction increases. This should 
force the Hausdorff dimension to decrease monotonously with decreasing E c on each ij. 

6.3 Multifractal spectrum and energy transport. 

The maps are not conformal. In this case, the multifractal spectrum is obtained from a 
(sub-additive) thermodynamic formalism [llj. The corresponding potential is a function of the 
singular values ai(u,k), i = 1 . . . N, k = 1 . . . oo, of the product matrix L^ k . . . . along a 
legal sequence u G £q. 

It is interesting to remark that these singular values are directly related to the eigenvalues of 
the tangent map -DFx of the relaxation map (0). On the other hand, one can shown that -DFx 
is the evolution operator of a Markov process that governs the propagation of energy in Zhang 
model [HIE]- The eigenmodes of -DFx for the singular diffusion (JIJ) are the analog of Fourier 
modes for normal diffusion. The eigenvalues of DFx define a hierarchy of characteristic times, 
where the effective energy transport has distinct characteristics: singular on short time scales, 
anomalous at intermediate scales, and normal on long time scales Consequently, there 

is a close connexion between the energy transport and the multifractal spectrum of /t^. Some 
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aspects of this connexion, involving Lyapunov exponents, have been discussed in j7J. Further 
developments are under investigation 0. 

6.4 Ledrappier- Young formula and critical exponent of avalanche 
sizes. 

The Ledrappier- Young formula jTSj relates the positive Lyapunov exponents, the corresponding 
partial Hausdorff dimensions and the Kolmogorov-Sinai entropy. Applied to the backward 
dynamics of Zhang model it writes: 

N 

- = log(A0 - log( J N ) (14) 

i=i 

where the criii) are the partial Hausdorff dimensions in the direction i. Let J/v(X) be the 
number of preimages of X then — j JAr(X)<i/i(X) is the averaged number of preimages. 
log(iV) — log( J/v) is therefore the backward entropy. In the case where the dynamics is invertible 
log( Jjv) = and the backward entropy is equal to the forward entropy. When the system is not 
invertible, one can make it invertible by coding the backward iteration tree in the same way 
as we did with the excitation sequences, hence introducing an additional variable on which the 
forward dynamics contracts. 

Except in some very restricted cases [5] log(iV)— log(Jjv) is proportional to log(iV). Since the 
partial Hausdorff dimensions obey u L (i) < 1 the equation (f*HJ) implies — J2i^L{i) > alog(iV). 
Then, from eq. (J12)) . we obtain an upper bound on the avalanche size: 

< s > L > —^--rlog(N) (15) 
I !og( e )l 

with strict equality iff all the partial dimensions are equal to 1. When fi^ is singular < s >l 
diverges therefore faster than logarithmically with N. This implies that the critical exponent 
of Pl{s) is such that r < 2. We get therefore a bound on the critical exponent of avalanche 
size from fractal considerations. 

7 Conclusion. 

We have shown that certain classes of models of SOC like the Zhang model fit naturally into 
a well known class of dynamical systems. Especially for the question of asymptotic energy 
distribution, observables distribution, ergodicity, this seems to be a proper point of view. Fur- 
thermore it exhibits close relationship between the probability of the size of avalanches and the 
fractality of the attractor. Further developments are under investigations. 
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